#include "../DWMRICore/Fiber.h"
#include "../DWMRICore/Point3I.h"
#include "../DWMRICore/Scalar.h"
#include "../DWMRICore/VectorField.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <vector>

using namespace std;

char *fiber_pathname = NULL;
char *dst_pathname = NULL;

void PrintUsage()
{
	printf("Usage: Fibers2VectorFields -fiber configure_file -dst dst_file\n");
}

int ParseArguments(const int argc, char *argv[])
{
	for (int i = 1; i < argc; ++i) {
		if (strcmp(argv[i], "-fiber") == 0) {
			fiber_pathname = argv[++i];
		} else if (strcmp(argv[i], "-dst") == 0) {
			dst_pathname = argv[++i];
		} else {
			return -1;
		}
	}

	if (fiber_pathname == NULL || dst_pathname == NULL)
		return -2;

	return 0;
}

int main(int argc, char *argv[])
{
	if (ParseArguments(argc, argv) != 0) {
		PrintUsage();
		return 0;
	}

	/* read the configure file */
	char **confPathname = NULL;
	char **vertPathname = NULL;
	int w, h, d, scale;
	float stepSize, maxLength;
	ReadConfigureFile(fiber_pathname, &confPathname, &vertPathname, w, h, d, scale, stepSize, maxLength);

	/* generate vector fields from fibers */
	vector<vector<CVector3F> > vectorFields;
	vectorFields.resize(w*h*d);
	int max_count = 0;
	for (int z = 0; z < d; ++z) {
		printf("\rGenerate vector fields from z=%d ... ", z);
		fflush(stdout);

		/* read fibers in current slice */
		CFiber **fibers;
		int *count;
		int cw, ch;
		ReadFibers(confPathname[z], vertPathname[z], &fibers, &count, cw, ch);

		for (int i = 0; i < cw * ch; ++i) {
			if (count[i] > max_count)
				max_count = count[i];

			for (int j = 0; j < count[i]; ++j) {
				vectorFields.at(z*w*h+i).push_back(fibers[i][j].m_seed.dir);
			}
		}

		for (int i = 0; i < cw * ch; ++i) {
			SafeDeleteArray(fibers[i]);
		}
		SafeDeleteArray(fibers);
		SafeDeleteArray(count);

		printf("Done.");
	}

	CVectorField *result = new CVectorField();
	result->CreateVectorField(max_count, w, h, d);
	for (int i = 0; i < w*h*d; ++i) {
		for (int j = 0; j < vectorFields.at(i).size(); ++j) {
			result->SetVector(i, j, vectorFields.at(i).at(j));
		}
	}
	result->SaveVectorFieldFile(dst_pathname);
	SafeDelete(result);


	for (int z = 0; z < d; ++z) {
		delete[] confPathname[z];
		delete[] vertPathname[z];
	}
	delete[] confPathname;
	delete[] vertPathname;

	return 0;
}


//void Fiber2VectorFields(CFiber *fiber, vector<vector<CVector3F> > &vectorFields, 
//	const int w, const int h, const int d, const int scale)
//{
//	vector<int> voxels;
//	voxels.clear();
//
//	int lastIndex = -1;
//	/* resample the fiber and add to the vector fields */
//	for (uint i = 1; i < fiber->m_fCount; ++i) {
//		CPoint3F sPt = fiber->m_pF[i-1];
//		CPoint3F ePt = fiber->m_pF[i];
//		CVector3F dir = ePt - sPt;
//		CVector3F ndir = dir;
//		ndir.Normalize();
//
//		for (float l = 0.0f; l < 1.0f; l += 0.1f) {
//			CPoint3F pt = sPt + dir * l;
//			int index = (int(pt.m_z) * scale * h + int(pt.m_y) * scale) * w + int(pt.m_x) * scale;
//			if (lastIndex == index)
//				continue;
//
//			int flag = -1;
//			for (uint k = 0; k < voxels.size(); ++k) {
//				if (voxels.at(k) == index) {
//					flag = 1;
//					break;
//				}
//			}
//			if (flag == -1) {
//				vectorFields.at(index).push_back(ndir);
//				lastIndex = index;
//			}
//		}
//	}
//	for (uint i = 1; i < fiber->m_bCount; ++i) {
//		CPoint3F sPt = fiber->m_pB[i-1];
//		CPoint3F ePt = fiber->m_pB[i];
//		CVector3F dir = ePt - sPt;
//		CVector3F ndir = dir;
//		ndir.Normalize();
//
//		for (float l = 0.0f; l < 1.0f; l += 0.1f) {
//			CPoint3F pt = sPt + dir * l;
//			int index = (int(pt.m_z) * scale * h + int(pt.m_y) * scale) * w + int(pt.m_x) * scale;
//			if (lastIndex == index)
//				continue;
//
//			int flag = -1;
//			for (uint k = 0; k < voxels.size(); ++k) {
//				if (voxels.at(k) == index) {
//					flag = 1;
//					break;
//				}
//			}
//			if (flag == -1) {
//				vectorFields.at(index).push_back(ndir);
//				lastIndex = index;
//			}
//		}
//	}
//}
//
//void MeanShift(vector<CVector3F> &in, vector<CVector3F> &out)
//{
//	if (in.size() < 2) {
//		for (uint i = 0; i < in.size(); ++i) {
//			out.push_back(in.at(i));
//		}
//		return;
//	}
//
//	vector<CVector3F> remain;
//	for (uint i = 0; i < in.size(); ++i) {
//		remain.push_back(in.at(i));
//	}
//	while (remain.size() > 0) {
//		CVector3F c = remain.at(0);
//		float delta = 0.0f;
//		do {
//			CVector3F new_c = CVector3F(0.0f, 0.0f, 0.0f);
//			for (uint i = 0; i < remain.size(); ++i) {
//				float t = InnerProduct(c, remain.at(i));
//				if (t > 0.9f) {
//					new_c = new_c + remain.at(i);
//				} else if (t < -0.9f) {
//					new_c = new_c - remain.at(i);
//				} else {
//					;
//				}
//			}
//			new_c.Normalize();
//			delta = 1.0f - InnerProduct(c, new_c);
//			c = new_c;
//		} while (delta > 0.001);
//		
//		out.push_back(c);
//
//		vector<CVector3F> temp;
//		for (uint i = 0; i < remain.size(); ++i) {
//			float t = InnerProduct(c, remain.at(i));
//			if (t > 0.9f || t < -0.9f) {
//				continue;
//			} else {
//				temp.push_back(remain.at(i));
//			}
//		}
//		remain.clear();
//		for (uint i = 0; i < temp.size(); ++i) {
//			remain.push_back(temp.at(i));
//		}
//	}
//}
//
//
//void SaveResults(vector<vector<CVector3F> > &result, const int w, const int h, const int d)
//{
//	FILE *fp = fopen(dst_pathname, "wb");
//	
//	fwrite(&w, sizeof(int), 1, fp);
//	fwrite(&h, sizeof(int), 1, fp);
//	fwrite(&d, sizeof(int), 1, fp);
//
//	for (int i = 0; i < w * h * d; ++i) {
//		int count = result.at(i).size();
//		fwrite(&count, sizeof(int), 1, fp);
//
//		for (int j = 0; j < count; ++j) {
//			CVector3F dir = result.at(i).at(j);
//			float v[3];
//			v[0] = dir.m_x;
//			v[1] = dir.m_y;
//			v[2] = dir.m_z;
//			fwrite(v, sizeof(float), 3, fp);
//		}
//	}
//
//	fclose(fp);
//}
//
//float ComputeVariance(vector<CVector3F> &in)
//{
//	for (uint i = 1; i < in.size(); ++i) {
//		float temp = InnerProduct(in.at(0), in.at(i));
//		if (temp < 0.0f) {
//			in.at(i).m_x = -in.at(i).m_x;
//			in.at(i).m_y = -in.at(i).m_y;
//			in.at(i).m_z = -in.at(i).m_z;
//		}
//	}
//
//	CVector3F mean(0.0f, 0.0f, 0.0f);
//	for (uint i = 0; i < in.size(); ++i) {
//		mean = mean + in.at(i);
//	}
//	mean.Normalize();
//
//	float cov = 0.0f;
//	for (uint i = 0; i < in.size(); ++i) {
//		cov += fabs(InnerProduct(in.at(i), mean));
//	}
//	cov = cov / (float)in.size();
//	return 1.0f/cov;
//}
//
//int main(int argc, char *argv[])
//{
//	if (ParseArguments(argc, argv) != 0) {
//		PrintUsage();
//		return 0;
//	}
//
//	/* read the configure file */
//	char **confPathname = NULL;
//	char **vertPathname = NULL;
//	int w, h, d, scale;
//	float stepSize, maxLength;
//	ReadConfigureFile(fiber_pathname, &confPathname, &vertPathname, w, h, d, scale, stepSize, maxLength);
//
//	/* generate vector fields from fibers */
//	vector<vector<CVector3F> > vectorFields;
//	vectorFields.resize(w*h*d);
//	for (int z = 0; z < d; ++z) {
//		printf("\rGenerate vector fields from z=%d ... ", z);
//		fflush(stdout);
//
//		/* read fibers in current slice */
//		CFiber *fibers;
//		int cw, ch;
//		ReadFibers(confPathname[z], vertPathname[z], &fibers, cw, ch);
//
//		for (int i = 0; i < cw * ch; ++i) {
//			Fiber2VectorFields(&(fibers[i]), vectorFields, w, h, d, scale);
//		}
//
//		SafeDeleteArray(fibers);
//
//		printf("Done.");
//	}
//
//	/* print some statistic information */
//	int average = 0;
//	int min = 1000;
//	int max = 0;
//	int max_index = 0;
//	for (int i = 0; i < w * h * d; ++i) {
//		average += vectorFields.at(i).size();
//		if (vectorFields.at(i).size() > max) {
//			max = vectorFields.at(i).size();
//			max_index = i;
//		}
//		if (vectorFields.at(i).size() < min)
//			min = vectorFields.at(i).size();
//	}
//	average = average / (w * h * d);
//	printf("\rmin = %d, max = %d, average = %d, max_index = %d\n", min, max, average, max_index);
//
//	/* compute the co-variance for each voxel */
//	CScalar *variance = new CScalar();
//	variance->CreateScalar(w, h, d);
//	for (int i = 0; i < w * h * d; ++i) {
//		int size = vectorFields.at(i).size();
//		if (size > 0) {
//			variance->SetScalar(i, ComputeVariance(vectorFields.at(i)));
//		} else {
//			variance->SetScalar(i, 0.0f);
//		}
//	}
//	variance->SaveScalarFile(dst_pathname);
//	delete variance;
//
//	///*  */
//	//vector<vector<CVector3F> > results;
//	//results.resize(w*h*d);
//	//for (int z = 0; z < d; ++z) {
//	//	for (int j = 0; j < w * h; ++j) {
//	//		int i = z * w * h + j;
//	//		MeanShift(vectorFields.at(i), results.at(i));
//	//	}
//	//	printf("\rz = %d", z);
//	//	fflush(stdout);
//	//}
//
//	//average = 0;
//	//min = 1000;
//	//max = 0;
//	//max_index = 0;
//	//for (int i = 0; i < w * h * d; ++i) {
//	//	average += results.at(i).size();
//	//	if (results.at(i).size() > max) {
//	//		max = results.at(i).size();
//	//		max_index = i;
//	//	}
//	//	if (results.at(i).size() < min)
//	//		min = results.at(i).size();
//	//}
//	//average = average / (w * h * d);
//	//printf("\nmin = %d, max = %d, average = %d, max_index = %d\n", min, max, average, max_index);
//
//	/* save */
//	//SaveResults(results, w, h, d);
//
//	for (int z = 0; z < d; ++z) {
//		delete[] confPathname[z];
//		delete[] vertPathname[z];
//	}
//	delete[] confPathname;
//	delete[] vertPathname;
//
//	return 0;
//}